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VAPORIZATION RESPONSE OF EVAPORATING DROPS 


WITH FINITE THERMAL CONDUCTIVITY 
by V. D. Agosta and S. S. Hammer 


SUMMARY 

The primary objective of the analysis was to obtain a numer- 
ical computing procedure for determining the vaporization response 
of droplets with finite thermal conductivity in an oscillating 
pressure and flow field. The governing equations for vaporization 
of liquid drops in a rocket combustor environment were taken from 
Refs, 1 and 2. Additional equations were employed to account for 
finite thermal conductivity of the liquid drop (Ref. 3) . The system 
of equations were solved utilizing a finite difference technique 
and a high speed digital computer. Oscillation in the rates of 
vaporization of an array of repetitively injected drops in the 
combustor were obtained from the summation of individual drop 
histories. A nonlinear in-phase frequency response of the entire 
vaporization process to pressure oscillations was calculated and a 
response factor, determined as defined by Equation 1 of 

Ref. 4. In addition, a nonlinear out-of-phase response factor, I^^ 
in-phase and out-of-phase harmonic response factors R^^, R 2 , 1^^, I 2 
and a Princeton type "n” and "t" were determined as described in 
Refs. 5 and 6. In general, it was found that the nonlinear in- 

phase response factor, R , was not very sensitive to variations 
3 . 

of up to 10 in the liquid thermal conductivity, for the frequency 
range of interest in combustion instability studies. 

The resulting data was correlated and is presented in graphical 
format. Qualitative agreement with the open literature is obtained 
in the behavior of the in-phase response factor. 
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INTRODUCTION 

Studies on nonlinear combustion instability have been per- 
formed (Refs. 2,4,7) which indicate that vaporizing drops are 
sensitive to frequency-dependent pressure oscillations. The sensi- 
tivity of the vaporization processes has been traced to thermal 
time lags, namely, that time delay between changes in the drop 
temperature and subsequent mass vaporization and changes in the 
drop environment. The thermal time lag model proposed in Ref. 7 
is extended to include the effects of drop thermal conductivity on 
drop surface temperaiiure and mass vaporization. 

Wave deformation effects on droplet vaporization as proposed 
in Ref. 5 were also considered in this investigation by varying 
the harmonic distortion in a pressure disturbance propagating in a 
liquid rocket combustion chamber. 

Nonlinear and harmonic in-phase and out-of-phase response 
factors which have been evolved from linear limit-cycle stability 
analyses have been adopted for use in this report. This pro- 
cedure is consistent with present usage and allows for the com- 
parison of data and results for similar parameters and factors. 


THEORY 


Drops of liquid oxygen are assumed to be vapdrizing in com- 
bustion gases, composed of stoichiometric reaction products with 
hydrogen in a cylindrical combustion chamber, with an established 
travelling- transverse acoustic mode. The instantaneous pressure 
and gas velocity fvinctions consist of the steady state and oscillating 
components attributed to the acoustic mode. They are expressed as 
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where 

a = speed of sound 

f = frequency of oscillations 

J = Bessel function 

m = mass of droplet 

m^= initial mass of droplet 

p = instantaneous pressure 

p^= mean chamber pressure 

APj^= peak-to~peak pressure amplitude (fundamental) 

AP 2 = peak-to-peak pressure amplitude (harmonic) 

R = radial location in chamber 

radius of chamber 

t = time 

T = temperature 

U = gas velocity 

U = final axial gas velocity 
^f 

0 = phase angle 
Y = isentropic exponent 

The wave is assumed to be adiabatic with chamber temperature 
and pressure related by 

■p .Y.~Ai 

=(:?-) '• ( 6 ) 

T p 

c -^c 
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The oscillations of pressure and velocity are superimposed 
on the mean level of parameters affecting drop evaporation and 
motion. The vaporization rate is controlled by heat and mass 
transfer to the surface of a drop with finite thermal conductivity. 
Drop motion is controlled by a momentum balance as a result of drag 
with the enveloping gas. Axial gas velocity is proportional to 
the fraction of drop mass vaporized, and it attains a final assumed 
velocity at complete evaporation. The complete drop history is 
defined by the following equations for weight evaporation rate 
heat transfer, acceleration in an axial direction and temperature 
distribution within the drop. 

The governing differential equation for temperature in a 
spherical liquid droplet is given by 



a(T 


rr 



0 < r < r 
— s 

t > 0 


(7) 


where a = liquid thermal diffusivity 

r = radial coordinate within the drop 

r = surface radius, 
s 

The subscripts r and t represent differentiation with respect 
to radius and time, respectively. For small droplets the assump- 
tion of spherical geometry is usually accepted. It is recognized 
that hwere drag exists, droplets deform; however, in this analysis 
the deformation is neglected. The initial condition assumes that 
the droplet is at uniform temperature T^, 

T(r,0) = T 0 < r * t=0 (8) 

V s 


The boundary conditions are 
T(0,t) is finite 


r=0 ; t>0 


and 
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where h = heat transfer coefficient 

Tg = drop surface temperature 

k = liquid thermal conductivity 

X = heat of vaporization 

Cr^ = specific heat of droplet vapor 
Pv 

w = mass evaporation rate 
= droplet surface area 

It is assumed that the film thickness surrounding the droplet is 
small compared to the droplet diameter. 

The mass evaporation is given by 


w = A^KgP 



r=r 

s 


( 11 ) 


where p^ = droplet vapor pressure 

K = mass transfer coefficient. 

g 

The heat and mass transfer coefficients h and K^, respectively, 
are obtained from the Nusselt correlations 
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T = 
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Re 


vapor-gas mixture thermal conductivity 
universal gas constant 
molecular weight of vapor 
arithmetic mean temperature, 

Prandtl number (c u/k) 

p^' 'm 


Tc+Ts 


= Schmidt number (jj/Dp) 
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drop velocity 


D = molecular diffusion coefficient 
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The droplet acceleration is determined by considering the 
momentum transfer, between the liquid drop and gases due to aero- 
dynamic drag 


^^d ^ 1 p hi ( AV) I AV 
dt 8 h) r^ 


and the drag coefficient 




27 


2p AVr 
"^m s 


-.84 


(16) 


where AV is the difference between the axial gas velocity and 
drop velocity. 

The droplet radius is, of course, a function of time and is 
related to the mass evaporation rate. 


r 

s 


w 




(17) 


where f = surface regression rate; time rate of change of drop 
s 

radius . 

The analysis of droplet evaporation in a gas stream, as form- 
ulated above, is developed into a computer program for the numerical 
solution of the time dependent evaporation rate, droplet radius and 
temperature distribution within the drop, A detailed discussion of 
the calculation procedure and a program listing are contained in 
Appendix B. Beginning with a specification of the initial conditions 
the droplet vaporization and surface regression rates are determined 
from Eqs, (11) and (17), The heat balance equation (10) at the drc?)- 
let surface is then used to determine the temperature gradient at 
the surface. By finite difference techniques the temperature gradient 
at a point adjacent to the surface, and the second derivative of 
temperature with respe<^to^r^^;t_the - surface-are determined, 
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Equation (7) is then used to obtain the variation of surface temp- 
erature with time. Interior point temperature calculations are 
performed by utilizing a finite difference scheme for the solution 
of Eq. (7) . At the end of the time interval new values are calcu- 
lated for surface radius, drop velocity, surface temperature and 
droplet mass by integrating ghe appropriate time derivatives over 
the interval. The gas pressvire, temperature and velocity are 
evaluated at the end of every interval from Eqs. (l)-(6). These 
thermodynamic and geometric properties then become the initial 
conditions for the next time interval . The procedure continues 
until the droplet mass is reduced to 10% of its initial value. 

A summation of single-drop histories is used to determine the 
time variations in vaporization rate of a one-dimensional array of 
repetitively injected drops. An arbitrary number of drops are 
injected into the chamber every cycle at times distributed evenly 
over the oscillation period. For the case of four drops injected 
per cycle these would appear in the chamber at intervals of one- 
quarter period. Vaporization histories vary among drops injected 
at different times during one pressure oscillation; however, drops 
injected at times one period apart experience identical, acoustic 
pressure and velocity fields and thus, have identical histories. 
Eventually the same number of drops are completely vaporized per 
cycle as are injected, and the number of drops in the array becomes 
constant over each full cycle. When this steady state condition 
is reached the fully developed array consists of a number of drops 
equal to the ratio of drop burning time to oscillation period times 
the number of drops injected per cycle. The droplets in the array 
have a decreasing mass down the length of the chamber and range in 
age from a new drop just injected to an old one almost completely 
vaporized. With the array fully developed, the mass vaporizing 
from the entire array of drops at any time is obtained from a 
summation of the vaporization of the individual drops that constitute 
the array. 
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since the history of each injected drop is calculated inde- 
pendently of the preceding or succeeding drops, the computation 
time is dependent on the number of drops injected per cycle. In 
order to simulate a continuous injection of drops, it would be 
necessary to analyze an unwieldy number of drops. Alternatively, 
it is possible to smooth out perturbations in the array vaporization 
history caused by the appearance and disappearance of a small 
finite number of drops by artificially inserting additional drops 
between those whose histories are calculated. The vaporization 
rates of the artificially injected drops are determined by inter- 
polation between the vaporization rates of the drops for which 
calculations are performed. Thus the instantaneous vaporization 
rate for the entire array is obtained from a summation of the 
calculated drop histories and interpolated artificial drop histories. 

The vaporization history of eight drops, injected at equal 
intervals during a pressure cycle is shown in Figure 1. The 
vaporization rate tends to be higher at both the maximum and min- 
imum pressure condition in the oscillation than at the mean pressure 
conditions. This is a result of lower total velocities at the mean 
pressure condition. For the case shown in Figure 1, the next drop 
injected (drop 9) would exhibit the identical behavior as drop 1, 

For a burning time of .001 seconds and a pressure oscillation 
period of .00033 seconds the fully developed array would consist 
of 30 drops with drops 1, 9, 17, 25 exhibiting identical behavior, 
drops 2, 10, 18, 26 exhibiting identical behavior, etc. The 
vaporization from the entire array is obtained by adding the evap- 
oration rates of all of the constituent drops at corresponding 
times in the pressure cycle. The results of this calculation, along 
with the pressure curve, is shown in Figure 2, . It is evident from 
this curve that the evaporation rate tends to be higher at times in 
the pressure cycle corresponding to maximum or minimum pressure 
with relatively lower evaporation rates occurring at the mean 
pressure . 
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The calculation shown in Figures 1 and 2 required 160 seconds 
of computing time on a CDC 6600 computer. 

In order to insure that certain computational procedures being 
used would yield results that were reproducible, the following tests 
were conducted: In the drop evaporation test runs, calculation was 

terminated after 90% of the drop mass had evaporated. In order to 
insure that this was not a premature cutoff point, one rvin was 
extended to the 97% mass evaporation point. In the former case, the 
response function was 0.548; in the latter, 0.545. The difference is 
negligible; thus a mass termination point of 90% was adopted. 

In the determination of the response factor, eight drops 

were introduced, evenly spaced during the pressure cycle. The 
cumulative mass evaporation from these drops was calculated and 
svibsequently the response factor. At low frequencies the drop array 
size was not statistically meaningful. The question posed was: 
are eight oxygen drops a statistical meaningful array at a frequency 
of 1000 Hertz? A quick answer was sought. In one run, 20 drops 
were used in a drop array and the resulting response factor thus 
determined did not vary significantly from the results for the 
eight-drop array. 

RESPONSE 

Imposing acoustic oscillations on the pressure and velocity 
field in a rocket combustion chamber causes perturbations in the 
vaporization rate of the array of droplets resident in the chamber. 

The oscillation in vaporization rate will exhibit harmonic components 
of the imposed frequency of acoustic oscillations. A series of 
response factors are calculated in this study in order to relate 
vaporization rate oscillations to pressure and velocity field 
oscillations. These response factors are the in-phase and out-of- 
phase components of the vaporization rate oscillations relative to 
the acoustic pressure oscillations, where both oscillations are 
given as fractional perturbations about the mean value of the variable. 
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can be extracted 


The nonlinear in-phase response factor, 
from the array vaporization rate and normalized by correlation 
procedure (Ref. 5) defined by 
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where W is the instantaneous evaporation rate from the entire 
array and W is the average value of W taken over a full cycle. 

The nonlinear out-of-phase response factor, is given by 
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Additional response factors were calculated as part of this 
study ; 
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Out-of-phase harmonic response factors 
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A Princeton type "n" and "t" defined in Ref. 6 are also cal- 


culated from the nonlinear response factors 
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RESULTS AND DISCUSSION 

A. Parametric Study of Response Functions 

A series of calculations were made, utilizing the previously 
described computer program, in order to determine the effect of a 
variety of boundary conditions on the magnitude of the response 
functions for evaporating liquid oxygen droplets. The thermodynamic 
properties used in the calculations were: 

-7 2 

Thermal diffusivity = 7.56x10 ft /sec 

3 

Liquid density =72 Ibm/ft 
Prandtl number = .712 
Gas viscosity = 4«2xl0”^ Ibm/ft-sec 
Vapor specific heat = .288 BTU/lbm°F 

— 5 o 

Liquid thermal conductivity = 2.21x10 BTU/ft-sec- F 
Gas thermal conductivity = 4.04x10"’^ BTU/ft-sec-°F 
Isentropic exponent = 1.135 
Initial drop temperature = 140°R 
Combustion gas temperature = 6500°R 

The vapor pressure and heat of vaporization of the liquid 
oxygen is evaluated as a function of the drop surface temperature by 

lA Q-3 1476 

2 16, OS'— m o cn 

p^(lbf/ft^) = e ^ 

UBTU/lbm) = 61.33 + .5916T - ,00248T^ 

The following parameters were varied over the range indicated: 

Chamber pressure: 100-600 psia 

Initial drop radius: 15-150 microns 

Final gas velocity: 400-2400 ft/sec 

Initial drop velocity: 50-200 ft/sec 

Frequency of oscillation: 200-30,000 cps 

Amplitude of fundamental ( Ap^^) ; .01p^-.8p^ 

Amplitude of harmonic (AP 2 ): 0-1.2 Apj^ 
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Table I gives the full range of tests conducted with boundary 
conditions and response functions. 

A cursory look at the values calculated for the in-phase non- 
linear response factor, indicates many values that are greater 

than those previously foxind by the methods of Refs, 2 and 4. The 
present study includes the effects of finite thermal conductivity 
which, as discussed previously, causes a significant decrease in the 
time required for the drop surface to reach an effective equilibrium 
temperature. In addition, the temperature of the drop center remains 
relatively unchanged over the entire lifetime of the drop with the 
temperature gradient in the drop concentrated at the drop outer 
surface. Both these factors make the drop more sensitive to fluctua- 
tions in ambient pressure and temperature with concomitant increased 
response factors. 

The major cause of larger response factors obtained in this 
study, however, appears to be due primarily to the effect of wave 
distortion, i.e,, inclusion of harmonics in the ambient pressure 
field. This causes more relative peaks and valleys in the impressed 
pressure oscillation with attendant increases in evaporation rate. 
Another significant contributing effect is the presence of the array 
of drops in a varying gas velocity field. 

Figures 3 and 4 indicate the quantitative effect of wave dis- 
tortion and amplitude on the nonlinear in-phase response factor, 

For the case of a constant ratio of the first harmonic to fxindamental 
wave pressure amplitude, AP2/^3^=0-8, Figure 3, the variation of 
R^^ is relatively independent of fvindamental amplitude provided 

> O.lP^. For the case where the fundamental amplitude is small, 
Ap^=0,01p^, the response factor increases significantly to a maximum 
value of R^^=1.68. This behavior agrees qualitatively with that 
shown in Figure 2 of NASA TN D-6287, Ref, 5, Figures 4a, 4b and 4c 
show the effect of variations in the relative magnitude of the harmonic 
amplitude. The maximum response curve occurs at a value of Ap2=0.8Apj^. 
This too is in agreement with the results found in Ref. 5. It is 
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observed from Figures 4a,b,c that the frequency at which the maximum 

value of R occurs decreases with increasing fundamental amplitude, 
nt 

Figure 4b contains the results of calculations made for zero 
harmonic content, AP2=0. As discussed above, the results obtained 
in Ref. 4, IJASA TN D-3749, also were for the case of no harmonic 
content in the pressure oscillation. The magnitude of the response 
factors calculated in this study for are lower than those 

obtained with non-zero harmonics and are consistent with those ob- 
tained in NASA TN D-3749, 


The final gas velocity greatly affects the value of R . This 

n ^ 

effect is observed in Figure 5 at all frequencies. A crossplot was 
made at 1200 Hz, Figure 6, and the final gas velocity extended on the 
low scale to 50 feet per second with drop velocities of 100 feet 
per second. It is seen that the value of the nonlinear in-phase 
response factor increases to 2.36 at a final gas velocity of 200 feet 
per second. It is also observed that the values of R^^ are about 
twice those observed in the previous Figures 3 and 4, In the regime 
values and behavior of the nonlinear response factor, 

R ,, agree both qualitatively and quantitatively with those shown in 
Figure 7, NASA TN D-6287, Ref, 5, 


An attempt was made to correlate all of the response factor data 
obtained utilizing the transformed frequency suggested in Refs. 2 and 
4, However, the buckshot nature of the resulting curves indicated 
that while qualitative agreement between the two studies is obtained, 
the previous correlations are unsuitable for the kind of model em- 
ployed in this study. In order to correlate the data it was necessary 
to use a transformed response factor together with a modified trans- 
formed frequency. The transformations were obtained by trial and 
error in an attempt to minimize the data scatter. It was found that 
rather than include the effects of wave amplitude and distortion in 
the transformed coordinates, the effect of these variables could best 
be seen by utilizing parametric curves. 
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The transformed frequency is defined as 



and the transformed response factor is 



Figure 7 shows the results of the correlation study with 
separate curves to depict the effect of wave distortion and amplitude. 
For fundamental amplitudes greater than -Ip^ the effect of funda- 
mental amplitude is not significant, but the effect of harmonic 
amplitude must be considered. Decreasing the fundamental amplitude 
to O.Olp^ causes significant increases in the response factor. 

The study was extended to include the correlation of the in- 
phase fundamental response factor and harmonic response factor R 2 / 
and these are shown in Figures 8 and 9. Both response factors are 
transformed similar to R^^ and correlated against the transformed 
frequency factor cited above. The fundamental response factor 
increases with the amplitude of the harmonic and is relatively in- 
sensitive to variations in the amplitude of the fundamental for 
Ap^^. Ip^. The magnitude of the harmonic response function increases 
with a decrease in the magnitude of the harmonic pressure oscillation 
component. Response factors, R 2 » on the order of 3.0 were calculated 
for harmonic amplitudes equal to 0.2 times the fundamental amplitude. 
The qualitative behavior, and the juxtaposition of these curves, 
agree with those shown in Figure 8, Ref, 5. The curves shown in 
Figures 7,8,9 represent the best fit through the available data. 
Deviations from the curve were generally less than 10% from the 
curve. 

Figures 10,11 and 12 represent the correlation of the results 
for the nonlinear out-of-phase response factor, the fundamental 
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out-of-phase response factore, and the harmonic out-of-phase 

response factor, I,. The I . curves were correlated with transforma- 
tions similar to those employed for the in-phase response factor. 

The shape of the curve is similar to a ciibic with a relative maximum 
point in the transformed frequency range of 300-400 cps and a rela- 
tive minimum at a frequency corresponding to the location of the 
maximum point on the curves. 

The results for the correlation of the out-of-phase fundamental 
response factor indicate that wave amplitude and distortion have no 
effect. In addition, excellent correlation was obtained by trans- 
forming by a multiple of (300/p ) and the frequency by a multiple 
3 /2 ^ 

of (r/50) . Thus the final gas velocity and initial drop velocity 

are not effective in varying the magnitude of 1^^. 

The correlation of the data for the out-of-phase harmonic 
response factor, required utilizing a transformed ordinate 



plotted against a transformed frequency 



Correlation also required the use of parametric curves to display 
the effect of wave distortion; i.e., the magnitude of I 2 decreases 
with increasing harmonic component amplitude. 

In Figure 13, a corelation of a Princeton type "n" was plotted 
vs. frequency factor, F. It is seen to remain essentially constant 
over a broad spectrum of frequency factor. The value of "n” increases 
with the amplitude of the harmonic pressure perturbation amplitude. 

Figure 14 is a plot of t vs. frequency of oscillation for a 
variety of drop radii, initial drop velocity, final gas velocity and 
amplitude o f pressure distu rbance. T he __var ia tion_of— the^value s- of 
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of the parameters cited above are those which normally occur in 
liquid rocket engine practice. The results indicate that it is possi- 
ble to correlate the values of t with the frequency of oscillation 
without any significant dependence on any of the other parameters. 

The correlation seems to indicate a relationship between t and f of 
the form, 

fT = constant. 

The value of t is determined from Eq. (26) . The right-hand side 

of the equation, i.e., R /I , has a value in the range 0.1 to 4.0 

n ^ ri ^ 

for the vast amjority of cases tested in this study. The correspond- 
ing values for 2nfT are in the range 250 to 350 degrees. These, 
when plotted for the specific case tested, yielded essentially a 
straight line on a log-log plot as is suggested in the above equation. 

The data from the above studies were regrouped for each of the 
seven response factors to show the effects of fundamental pressure 
perturbation amplitude ( Ap^) , harmonic pressure perturbation amplitude 
(AP 2 )» injection velocity (Vj^) , final gas velocity (u^) , droplet 
radius (r) and chamber pressure (P^) . The following table shows a 
summary of the relative effects of the various parameters on the 
different response factors; 



APi 

CN 

% 


^f 

R 

P 

c 

R 

weak 

moderate 

weak 

strong 

strong 

weak 

weak 

strong 

weak 

strong 

strong 

moderate 

«2 

weak 

strong 

weak 

strong 

strong 

weak 


weak 

moderate 

weak 

moderate 

moderate 

weak 

^1 

weak 

weak 

weak 

weak 

moderate 

weak 

^2 

weak 

moderate 

moderate 

moderate 

moderate 

weak 

n 

weak 

strong 

moderate 

strong 

moderate 

weak 

T 

independent of all parameters 
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B. Effect of Thermal Conductivity 


In general, it was fotind that the nonlinear in-phase response 

3 

factor, R , was not very sensitive to variations of up to 10 
n t 

in the liquid thermal conductivity, for the frequency range of 
interest in combustion instability studies. These results are 
explained below. 

The calculation for the R requires the summation of the 

nt ^ 

instantaneous evaporation rate of the complete array of drops in the 
chamber over an interval of time. Thus, any particular local event 
resulting from the interaction of the pressure wave and one drop is 
masked by the summation of all the other events occurring at the 
same time. In other words, the random event associated with each 
drop is masked by the stoichastic behavior of the drop array. \ 


A second reason for the insensitivity of R^^ to the thermal 
conductivity variation is that an open loop analysis is assumed. 

There is no feedback to the wave behavior from the droplet evapora- 
tion, This feedback can be significant in deforming the wave shape, 
and consequently the response factor, The wave shape can be 

deformed not only by coupling the relaxation time to the evaporation 
but also by the amount of vaporization occurring from the array of 
drops during the coupling. These two situations are discussed below 
for a particular case. 

The drop is introduced at a uniform temperature of 140°R. The 
drop surface heats up to about 240°R and remains at that temperature 
during its lifetime. For a 50p radius oxygen drop, and normal thermal 
conductivity, it takes 1,99x10 sec. for the drop surface to reach 
238°R, For the drop with high thermal conductivity, i.e., 10 ^ Btu 
ft ^°R ^sec it takes 1.87x10 ^ sec. for the drop surface to reach 
239°R, It is noticed that the time lag differs by about two orders 
of magnitude. In a closed loop analysis, this difference in time lag 
could mean the difference in the effective coupling of the drop 
evaporation to drive the passing pressure wave. 



For the case where the normal thermal conductivity is employed, 
the drop center does not appreciably heat up during its lifetime. 

For example, after 90% of the drop mass has evaporated, the temp- 
erature of the drop center is still 144°R. As found previously 
(Ref. 3), the temperature gradient is concentrated at the drop 
outer surface. For the case where the large thermal conductivity 
is used, it is found that the drop remains uniform in temperature. 

The uniform temperature drop has a much larger thermal inertia 
than the normal drop; thus the dynamic behavior of these drops would 
be very different, not only in the lag times, but also in the mass 
evaporation response to a disturbance. For example, for a passing 
compression wave, the normal drop would obtain a much greater surface 
temperature than the uniform temperature drop for the same amount 
of heat transfer to the droplet. If one includes the nonlinearity 
of the Clausius-Clapeyron equation, then the additional mass evap- 
oration in the former case would be much greater, with increased 
concomitant response factor. 

SUMMARY OF RESULTS 

The objective of this investigation was to analytically determine 
stability parameters relating drop vaporization response rates in a 
liquid rocket combustor to nonlinear-high amplitude pressure oscilla- 
tions for drops vaporizing with finite liquid thermal conductivity. 

The results of the program can be summarized as follows: 

1. A computer program was developed for determining the vapor- 
ization response of droplets with finite thermal conductivity to high 
amplitude distorted pressure oscillations. 

2. Drop vaporization responses for an array of drops traveling 
through a rocket combustor are significantly different than the 
response of a single drop stationary in a flow field. 

3. The vaporization rate responses were very dependent on drop 
size, final gas velocity in the combustor and wave distortion. The 
responses were moderately dependent on liquid droplet velocity and 
weakly dependent on chamber pressure and wave amplitude. 
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4. A general correlation scheme was developed to relate 
various response numbers that can be used in stability analyses to 
combustor operating conditions. 

5. For the Princeton type time lag response, the time lag was 
found to be inversely proportional to the frequency, resulting in a 
phase lag of 250 to 350 degrees. 
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TABLE I 


*PARAMETRIC STUDY OF RESPONSE FUNCTIONS 


Unless otherwise noted, the bo\indary conditions used to calculate 
the response functions are: 

Chamber pressure p^ = 300 psia 
Drop radius = 50 microns 
Drop initial temperature = 140°R 
Drop initial velocity = 100 ft/sec 
Final gas velocity = 1200 ft/sec 

Amplitude of fundamental pressiire perturbation; Ap^=.2 p^ 
Amplitude of harmonic pressure perturbation; AP 2=.8 APj^ 


f 


^1 

^2 


^1 

^2 

n 

TXIO^ 

150 

.289 

.244 

.359 

.345 

.252 

.491 

.404 

5.7 

200 

.351 

.276 

.468 

.442 

.389 

.524 

.519 

4.3 

400 

.581 

.526 

.666 

.448 

.606 

.20 

.541 

2.019 

800 

.810 

.910 

.654 

.185 

.265 

.0603 

,538 

.873 

1500 

.986 

1.112 

.789 

.232 

.216 

.025 

.65 

.46 

3000 

.941 

1.04 

.781 

.255 

.193 

.352 

.682 

.24 

6000 

.831 

.987 

.628 

.44 

.355 

.583 

.63 

.125 

12000 

.549 

.8 

.255 

.684 

.609 

.802 

.73 

.071 

30000 

.086 

.340 

-.309 

.765 

.808 

.967 

.4 

.032 

APj^^=.l 









600 

.674 

.720 

.602 

.425 

.558 

.216 

.555 

1.3 

1200 

.911 

.854 

1.0 

.131 

.189 

.0415 

.606 

.56 

2400 

1.026 

1.032 

.823 

.241 

.162 

.411 

.632 

.316 

5000 

.844 

.979 

.635 

.412 

.329 

.542 

.627 

.152 

10000 

.623 

.843 

.275 

.630 

.547 

.760 

.729 

.084 


21 



f 

R . 
n-C 

^1 

^2 


"l 

^2 

n 

TxlO^ 

APj ^=.4 

P 

c 








250 

.431 

.316 

.612 

.438 

.466 

.394 

.525 

3.35 

500 

.724 

.783 

.631 

.331 

.503 

.061 

.527 

1.507 

1200 

.891 

.998 

.725 

.160 

.214 

.071 

.590 

.57 

3000 

.920 

1.012 

,757 

.240 

.182 

.338 

.615 

.236 

8000 

.742 

.916 

.470 

.518 

.423 

.666 

.647 

.099 

I-* 

II 

• 

CD 

P 

c 








400 

.683 

.672 

.7 

.322 

.485 

.067 

.501 

1.89 

800 

.948 

1.091 

.724 

-.001 

.003 

-.0072 

.674 

.84 

1600 

.941 

.979 

.882 

.115 

.0625 

.197 

.628 

.41 

2500 

.918 

.995 

.797 

.196 

.152 

.265 

.609 

.277 

6000 

.814 

.938 

.621 

.387 

.299 

.526 

.599 

.126 

o 

II 

CN 

APj ^=. 

2 P 

c 







800 

.566 

.566 

— 

.075 

.075 

__ 

.3 

.67 

1500 

.692 

.692 

- 

.096 

.096 

- 

.305 

.331 

3000 

.488 

.488 

- 

.176 

.176 

- 

.273 

.202 

6000 

.43 

.43 

- 

.32 

.32 

- 

.335 

.117 

12000 

.244 

.244 

- 

.54 

.54 

- 

.673 

.072 

AP 2=.2 









800 

.791 

.656 

3.16 

.102 

.137 

.322 

.418 

.68 

1600 

.874 

.722 

3.61 

.185 

.176 

.399 

.470 

.36 

3000 

.72 

.635 

3.16 

.188 

.209 

.580 

.400 

.196 

6000 

.655 

. 566 

2.95 

.346 

.315 

.792 

.430 

.111 

12000 

.455 

.432 

2.5 

.530 

.525 

1.11 

.550 

.067 
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f 

R 

R, 

R^ 

I . 

I, 


n 

TxlO 


n-t. 

1 

2 

nZ 

1 

2 



4P2=1.2 

APi 








400 

.548 

.618 

.499 

.334 

.716 

.068 

.427 

2.03 

800 

»689 

1.04 

.445 

.130 

.295 

.015 

.500 

.448 

1500 

.87 

1.30 

.574 

.250 

.193 

.264 

.557 

.27 

3000 

,831 

1.33 

.587 

.362 

.302 

.403 

.538 

.254 

6000 

.698 

1.23 

.408 

.461 

.371 

.509 

.583 

.133 

12000 

.380 

.984 

.098 

.710 

.627 

.755 

.941 

.076 

to 

II 

• 

to 

APi 

AP^=.l 

P 








600 

.642 

.551 

2.91 

.411 

.352 

.191 

.470 

1.155 

1500 

.873 

.743 

4.135 

.280 

.291 

-.005 

.500 

.404 

4000 

.821 

.725 

3.195 

.151 

.156 

.472 

.441 

.140 

10000 

.519 

.441 

2.475 

.536 

.516 

.665 

.557 

.077 

^2=1 ‘2 

APi 

APj^=.l 

P 

c 






600 

.580 

.788 

.435 

.302, 

.630 . 

.075 

.422 

1.327 

1200 

.773 

1.067 

,568 

.205 

.300 

.139 

.644 

.325 

2500 

.956 

1.532 

.62 

.33 

.060 

.437 

.719 

.317 

5000 

.749 

1.23 

.410 

.388 

.131 

.494 

.544 

.159 

10000 

.475 

1.051 

-.091 

.664 

.696 

.732 

.780 

.088 

AP2=1 . 2 

APi 

AP^=.8 

P 

c 






400 

.631 

.774 

,531 

.034 

.538 

.029 

.418 

1.93 

1000 

.922 

1.146 

,513 

.110 

-.2 

.128 

.486 

.557 

2500 

.831 

1.098 

.575 

.174 

-.214 

.216 

.487 

.255 

6000 

.720 

1.12 

.539 

.415 

.167 

.496 

.512 

.134 

AP2=.2 

APi 

H 

II 

: *CD 

P 

c 






400 

.558 

,498 

2.06 

.385 

.378 

.557 

.428 

1.76 

800 

.798 

.722 

2.69 

-.041 

,0752 

.229 

.420 

.603 

1600 

.774 

.688 

2.99 

+ .032 

-.0429 

.100 

.368 

.322 


23 



f 

R 



I . 


I-. 

n 

TXlO 


nl 

1 

2 

nt 

1 

2 



2500 

.784 

.671 

2.52 

.098 

-.0166 

.119 

.379 

.217 

6000 

.688 

.61 

2.27 

.2926 

.0934 

.632 

.360 

.106 

to 

II 

• 

CD 

APi 

AP ^=.2 

P U = 

c g 

=400 





300 

1.43 

1.636 

1.112 

.848 

1.130 

.403 

1.147 

2.6 

600 

1.950 

2.22 

1.510 

.3177 

.3677 

.239 

1.290 

1.13 

1000 

2.082 

2.50 

1.429 

.090 

.344 

.463 

1.380 

.68 

2000 

1.720 

2.12 

1.340 

.182 

.0805 

.366 

1.166 

.329 

4000 

1.645 

1.941 

1.181 

.353 

.262 

.496 

1.095 

.1755 

9000 

1.390 

1.755 

.821 

.616 

.502 

.786 

1.007 

.0834 

U =2400 
9 








400 

.414 

.322 

.557 

.416 

.466 

.326 

.476 

2.09 

800 

.608 

.636 

.563 

.260 

.391 

. 056 

.435 

.935 

2000 

.699 

.743 

.630 

.086 

.1 

.062 

.466 

.332 

5000 

.693 

.776 

.562 

.383 

.330 

.466 

.537 

.152 

12000 

.425 

.603 

.147 

.595 

.518 

.715 

.719 

.076 

U =400 
9 

II 

U1 

o 







300 

1.11 

1.185 

.994 

.715 

.977 

.304 

.926 

2.62 

600 

1.533 

1.750 

1.196 

.162 

.262 

.0056 

1.028 

1.098 

1500 

1.724 

1.933 

1.397 

.147 

.126 

.180 

1.164 

.434 

4000 

1.551 

1.8 

1.164 

.231 

.178 

.315 

1.029 

.168 

9000 

1.298 

1.63 

.766 

.587 

.487 

.743 

.941 

.0836 

U =2400 V„= 

9 D 

50 







400 

.398 

.313 

.531 

.430 

.440 

.347 

.465 

2.09 

800 

.599 

.623 

.561 

.244 

.390 

.0179 

.424 

.929 

2000 

.679 

.663 

.575 

.035 

.066 

-.0133 

.431 

.319 

5000 

.652 

.724 

.539 

. 359 " 

.318 

.422 

.505 

.154 

12000 

.406 

.578 

.138 

.591 

.518 

.703 

.724 

.073 
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f 

R . 
n-t 

^1 

^2 


^1 

^2 

n 

TXlO^ 

U =400 

. g _ 

1 Vj ^=200 







600 

2.415 

2.876 

1.694 

.254 

.246 

.266 

1.619 

1.097 

1000 

2.00 

2.59 

1.064 

.145 

.075 

.255 

1.355 

.648 

2000 

1,805 

2.036 

1.438 

-.026 

-.047 

.0057 

1.275 

.191 

4000 

1.707 

1.96 

1.35 

.135 

.06 

.253 

1.155 

.168 

9000 

1.5 

1.86 

.937 

.515 

.382 

.723 

1.03 

.0809 

U «2400 V ^= 

^ . D 

200 







800 

.716 

.868 

.478 

.374 

.402 

.332 

.544 

.958 

2000 

.779 

.939 

.528 

.237 

.101 

.451 

.527 

.359 

5000 

.700 

.758 

.610 

.339 

.2811 

.430 

.518 

.152 

12000 

.442 

.633 

.148 

.601 

.523 

.721 

.718 

.072 

U =1200 V = 
D 

200 







200 

.441 

.331 

.612 

.499 

.459 

.564 

.578 

4.25 

800 

1.056 

1.087 

1,008 

.209 

.308 

.0557 

.7 

.862 

1500 

1.22 

1.44 

.867 

.1482 

.519 

.425 

.858 

.494 

3000 

1.019 

1.2 

.735 

.288 

.131 

.535 

.685 

.237 

6000 

.935 

1.101 

.675 

.425 

.334 

,567 

.679 

.125 

12000 

.586 

.85 

,173 

.695 

.621 

.816 

.810 

.071 

U =1200 V= 

■ g ■ D 

50 







200 

.321 

.253 

.426 

.416 

.354 

.511 

.494 

4.32 

800 

.733 

.82 

.596 

.161 

.223 

.0646 

.486 

.87 

1500 

.825 

.889 

.724 

.084 

.113 

,037 

.554 

.44 

3000 

.923 

1.004 

.795 

.323 

.305 

.350 

.635 

.243 

6000 

.79 

.928 

.574 

-452 

.384 

.557 

.622 

.129 

12000 

.508 

.738 

.149 

.667 

.597 

.777 

.791 

.072 
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f 

R 

R-. 


I . 

I , 

It 

n 

TXlO 


ni 

1 

2 

ni 

1 

2 



r = 15 n 

U =1200 V = 
D 

100 






2000 

.431 

.340 

.573 

.473 

.507 

.428 

.551 

.423 

4000 

.660 

.713 

.577 

.374 

.409 

.162 

.518 

.193 

7000 

.792 

.822 

.745 

.197 

.225 

.153 

.528 

.1 

15000 

.811 

.914 

.627 

.474 

.18 

.538 

.641 

.05 

30000 

.578 

.781 

.261 

.600 

.325 

.719 

.695 

.028 

r = 150 U 








150 

.888 

.990 

.730 

.168 

.259 

.0266 

.589 

4.58 

300 

1.111 

1.25 

.913 

.299 

.17 

.228 

.747 

2.31 

600 

1.052 

1.217 

.793 

.462 

.297 

.563 

.758 

1.25 

1200 

.822 

1.029 

.497 

.575 

.501 

.691 

.718 

.664 

3000 

.261 

.565 

-.233 

.740 

.718 

.665 

1.342 

.311 

r = 50 l ^ 

p =100 







500 

c 

.712 

.797 

.580 

.195 

.258 

.096 

.477 

1.42 

1000 

.853 

.967 

.676 

.162 

.172 

.147 

.565 

.687 

2000 

.870 

.953 

.741 

.191 

.165 

.232 

.577 

.348 

4000 

.783 

.897 

.603 

.281 

.211 

.390 

.541 

.183 

8000 

.629 

.799 

.362 

.523 

.429 

.669 

.620 

.102 

p =600 









500 

.568 

.498 

.677 

.501 

.637 

.296 

.587 

1.65 

1000 

.847 

.959 

.672 

.255 

.348 

.109 

.572 

.718 

2000 

1.05 

1.167 

.866 

.3 

.285 

.325 

.706 

.357 

4000 

.974 

1.105 

.769 

.343 

.262 

.469 

.671 

.183 

8000 

.803 

1.00 

.493 

.567 

.475 

.711 

.706 

.099 

16000 

.401 

.711 

-.083 

.792 

.760 

.842 

1.119 

.056 
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Symbols 


s 

D 

f 

h 

k 

M 

m 

Pc 

Pv 

AP]_ 

AP2 

r 

Re 

t 

T 

T 

•^s 

U 

Yd 

w 

a 

Y 

P 

X 


speed of sovind 

droplet surface area 

drag coefficient 

droplet vapor specific heat 

molecular diffusion coefficient 

frequency 

film heat transfer coefficient 

liquid thermal conductivity 

vapor-gas mixture thermal conductivity 

mass transfer coefficient 

molecular weight 

droplet mass 

mean chamber pressure 

vapor pressure 

amplitude of fundamental pressure perturbation 
amplitude of harmonic pressure perturbation 
droplet radius 

regression rate of droplet surface 

Reynolds number 

time 

chamber temperature 
droplet surface temperature 
gas velocity 
drop velocity 
mass evaporation rate 
thermal diffusivity 
isentropic exponent 
density 

heat of vaporization 
viscosity 
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APPENDIX B 


calculation Procedure and Program Listing 

The procedure used in calculating the droplet evaporation 
histories, temperature distribution within the drop and response 
factors is illustrated by the flow diagram. 

Step 1, Load into the machine the following boundary condi- 

tions, initial conditions and computational parameters. 

Card 1, (14 15) NRUN - Number of test case 

NMONTH - Month 

MDAY - Day 

MYEAR - Year 

JA - Number of calculation steps between 

output for each drop 

NA - Nximber of mesh points minus one within 

drop 

NP - Number of drops injected per cycle 

NY - Number of summation histories per period 

NART - Number of artificial drops insisted 

between each of tne NP calculated drops 

Card 2, (8E 10.4) S - Initial drop radius (ft.) 

2 

PO - Mean chamber pressure (Ibf/ft ) 

VGAF - Final gas velocity (ft/sec) 

VDI - Initial drop velocity (ft/sec) 

DPC - Ratio of peak-to-peak f\indamental pressure 
oscillation to mean chamber pressure 
DPCI - Ratio of harmonic pressure oscillation 
to fundamental oscillations 
OMEGA - Frequency of pressure oscillations (cps) 

A - Stretching parameter (1.3) 
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Card 3. (8E 10.4) STAB 

TFO 


- Initial time increment (sec.) 

- Mean combusticai gas temperatvire (°R) 

TO - Initial droplet temperatvire (°R) 

THETA - Phase angle of pressure oscillation 

1 841R 

AJq - Zero order Bessel function Jq (■ ’ ) 

1 84lR 

AJj^ - First order Bessel fvinction Jj^ 

1 841R 

AJ 2 - Second order Bessel function J 2 (■ ' *^^ "'“ "0 

RWR - Ratio of radial location in chamber to 
wall radius 

Card 4. (8E 10.4) GAMMA - Ratio of specific heats for combustion 

gases 

2 

ALPHA - Liquid thermal diffusivity (ft /sec) 

RHOL - Liquid density (Ibm/ft^) 

RR - Universal gas constant 

2 

PCA - Critical pressure - droplet (Ibf/ft ) 

PCB - Critical pressure - combustion gases 

(Ibf/ft^) 

TCA - Critical temperature - droplet (°R) 

TCB - Critical temperature - combustion gases (°R) 

Card 5. (8E 10.4) EMA - Molecular weight - droplet 

EMB - Molecular weight - gases 

AA - Constant in diffusivity equation - 

2.33 X 10 "^ 

BB - Constant in diffusivity equation - 1.823 

EMV - Molecular weight of drop vapor 

PR - Prandtl number of gases 

AKB - Combustion gas thermal conductivity 

(BTU/ft-sec-°R) 

CPV - Drop vapor specific heat 

card 6. (8E 10.4) VIS - Combustion gas viscosity (Ibm/ft-sec) 

AK - Drop thermal conductivity (BTU/ft-sec-°R) 
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step 2 


wh 


Step 3 


Step 4 


Write out the value of all input variables and cal- 
culate all parameters that are constant throughout calcula- 
tions. Initialize droplet and qas parameters to begin 
calculation of N drop history. Redefine variables as 
follows: 

U = Tr (Al) 

so that Eq. (7) becomes 

(A2) 


The space variable within the drop is also redefined: 

a = r/s (A3) 

where s is the drop surface radius as a function of time. 
The space variable is further redefined as: 

X = Aa(z+l/(A-z)) - a (A4) 

a = (A-1 )/(a2+1-A) (A5) 


and A is a stretching parameter which concentrates mesh 
points near the surface of the drop. The heat conduction 
equation (A2) can now be written: 


oS . zY 5 „2 

V ‘-f- ^ > V ^ “xx 


(A6) 


where Y 


dx 
dz • 


The interior point calculations are performed 
(temperature distribution within the drop) by utilizing 
a finite difference scheme for the solution of the trans- 
formed heat conduction equation based upon the current 
values of drop surface temperature and heat transfer rate 
to the drop. 

Droplet evaporation rates, changes in combustion 
gas properties, droplet velocity and net heat transfer 
to the drop are calculated by solving Eqs. (10) through 
(17) . The evaporation rate is stored for future summation 
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in the array history, intermediate results are printed 
at the end of every JA time increments. 

Step 5. If the drop is not 90% evaporated, time and other 

thermodynamic and dynamic variables are incremented and 
control is transferred to Step 3. 

Step 6. If the drop has reduced to 10% of its initial mass, 

the evaporation rate is transferred into the array summa- 
tion matrix at a predetermined number of points during 
the period (NY) . 

Step 7. An interpolation of the evaporate rate between the 

th til 

N and (N-1)^ drop is made for each of the artificial 
drops injected between two drops for which calculations 
are performed. These evaporation rates are added to the 
array matrix. 

Step 8. Test for end of evaporation calculations. 

Step 9. Calculation of response factors by solving Eqs. (18) 

through (26) . 

Step 10. Print results. , 

The machine computation time is primarily a function of the 

number of drops injected per cycle. For the case described in 

Figs. 1 and 2, i.e., eight drops injected per cycle, a total of 

160 seconds of CDC 6600 computer time was required. 
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PROGRAM DROPS(INPUT.OUTPUT»TAPE5 = INPUT»TAPF6 = OUTPUT) 

C 

C PROGRAM DROPS WITH RNL ♦ INL « R 1 » R2» II * 1 2 *F ♦ G ♦ AND N AND TAU AUG. 7 

C 


000003 


COMMON 
IMl .SIT 

000003 


COMMON 

000003 

C 

COMMON 

1R(64) 

000003 


PV(W) : 

000016 


SL(W) : 

000025 


D(W.O) 


WWW (3000) «TTT (3000) «W (400) . 1 NDFX , KK2 * VEL ♦ VO . AKG . DO» PPV ♦ H 


T 


V>AVW(J0«40U) 

X(64) .2(64) iA8 (64 )•« AC (64) . AO (64 ) . UO (64 ) , U (64 ) . UM ( 64 ) . T (64 ) 


2980. 957qfl67*EXP( a. 928-1476. 5/ (W-3. 568 ) ) 
61 .33+.5916*W-.00248*W**2 
= 0D1/W*AA* (0/D02)**BB 


C 

C NA=NUMBER OF MESH POINTS MINUS ONE 

C JA=NUMBER OF STEPS BETWEEN OUTPUTS FOR THE FIRST DROP 

C NP=NUMBER OF DROPS INJECTED PER CYCLE 

C. NY=NUMBER OF SUMMATION HISTORIES PER PERIOD (400 MAXIMUM) 

C A=STRETCHING PARAMETER 

C 


000036 

000064 

000110 

000134 

000160 

000204 

000214 

000221 

000241 

000244 

000247 

000253 

000254 

000256 

000260 

000262 

000264 

000266 

000271 

000303 

000335 

000357 

000431 

000435 

000447 


READ (5.250) NRUN. MONTH «MDA Y, MYE AR. JA.NA.NP. NY. NART 
READ (5. 260) S.PO . VGAF . VDl .DPC.OPCl .OMEGA. A 
READ (5. 260) ST AB , TFO ♦ TO .THETA. A JO . A J1 ♦ A J2 . RWR 
READ (5.260) GAMMA. ALPHA. RHOL.RR.PCA.PCB.TCA.TCB 
READ (5.260) EMA.EMB.AA.BB.EMV.PR.AKB.CPV 
READ(5.260) VIS.AK 
DD2 = SORT (TCA«TCB) 

ODl = (PCA*PCB)*».333333/SQRT(EMA*EMB/(EMA-.EMB) )*(TCA«TCB)«*.4166^ 
17 

TAV = (TF0+T0)*.5 
ROAV = P0*EMB/ (RR*TAV) 

SCl-t = V 1 3/KUfl V/U t PU . I A V ) 

SIT = S 
NC = NA+1 
TAU = 1. /OMEGA 
NX = NY-1 

DTST = TAU/FLOAT (NX) 

DTBAR = TAU/FLOAT (NP) 

WRITE (6.320) 

WRITE(6.330)NP.NART.NY 

WRITE(6.270)NRUN.MONTH.MDAY.MYEAR. JA.PO.TFO.S.TO.VGAF.VDI 
WRITE (6.340) OMEG A . THET A. DPC. A JO . A J1 . A J2 . RWR 

WRITE (6.280) ALPHA.RHOL.RR.PCA.PCB.TCA.TCB.EMA.EMB.EMV. AA.BB.SCH 
IPR.AKB.CPV.AK. GAMMA. VIS 
WRITE (6.290) 

WRITE (6.300) NA.A.STAB 
WRITE(6.410)DPC1 


000455 

000457 

000460 

000464 

000466 


C 

C 

C 

C 

C 


CONVERT DPC TO 1/2 PEAK TO PEAK PERTURBATION 

CONVERT DPCl TO RATIO OF DPI TO PO 

DPC=0PC/2. 

DPC1=DPC*DPC1 
PI = 4.*ATAN(1.) 

DUM8=PI*4.*RH0L 
PI2 = 2.*PI 
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000467 

000471 

C0047.1 

000475 

000477 

000500 

000504 

000510 

000513 

000515 

000516 

000520 

000521 

000526 

000527 

000530 

000533 

000534 

000536 

000537 

C 

C 

C 

000547 

000550 

000552 

000553- 

000554 

000555 

000556 

000557 

000560 

000561 

000562 

000566 

000567 

000570 

000571 

000613 

000620 

000623 

000626 

000632 

000636 

000637 

000642 

000656 

000660 

000665 

000666 

000667 

000672 

000675 

000700 

000701 

000704 

000707 


PID2 = Pl/2. 

P0ME5A = PI2*0MFGA 
GA = (GAMMA-1 .) /GAMMA 
VD = VDI 
TI = 0. 

CCl=nPC« (A JO- A J21*. 430/GAMMA 

CC2=RWR*AJ1 *.467*DPC/GAMMA 

0UM4 = EMV/2./RR 

DX = l./FLOAT(NAl 

ODX = l./OX 

OnX2 = l./DX/OX 

TW0DX=2.*nX 

SIG = (A-1.)/(A«(A-1.)*1.) 

Al = A*SIG 
A2 = 2.*A1 
A3 = SIG*(A**2*1.) 

A4 = A2«*2/SIG 
DO 1 1 = 100 
DO 1 J=l*400 
1 SAVWn«J)=0. 

IMITIALIZATION OF VALUES FOR THE (NN)TH DROP 

DO 220 NN = 1,NP 
S = SIT 
S2 = S**2 
SN = S 
WDOT = 0. 

STS = 0. 

ST = 0. 

K = 0 
J = 0 
OT = STAB 

TIME = DT8AR*(NN-1) 

KKl = 0 
KK2 = 0 
VO = VDI 

P=P0* (1 .♦.B59eDPC*AJl*C0S(P0MEGA*TIME-THETA) ♦.B59*DPC1*AJ1*C0S(2. 
IPOMEGAoTIMEM 
TF = TF0» (P/P0)«*GA 
TAV = (TFO+TO)*.5 
ROAV = P0*EM8/(RR«*TAV) 

SCH = VIS/ROAV/D (P-TAV) 

0UM2 = ,6*PR«*. 333333 
DO 10 N = 1,NC 
X (N) = DX*FLOAT (N-1 ) 

Z(N) = (A3+X(N)-S0RT((A3*X(N))«*2-A4*X(N)))/A2 
R(N) = Z(N)*S 
F = X(N)/A1+1,/A-Z(N) 

F? = F«F 
F3 = F«F2 
Y = A1*(1.*F2) 

AB(N) = Z(N)«Y 
AC(N) = A2*F3«ALPHA 
T (N» = TO 

IF (N.EQ.A/A) YNA = Y 
AD(N) = Y»*2*ALPHA 
U(N) = R(N)«T0 
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00071 1 
0007 j»3 


000715 

000726 

000730 

000731 

000733 

000735 

000737 

000741 

000742 

000746 

000751 

000751 

000754 

000756 


000757 

000760 

000762 

000764 

000766 

000771 

000776 

001002 

001004 

001007 

001011 

001015 

001016 

001024 


001027 

001030 

001031 

001032 

001033 

001035 

001035 

001041 

001050 

001053 

001054 

001057 

001060 

001061 

001064 

001106 

001113 

001116 


10 CONTTNUF 
TINT = TIMF 
C 

C TIME DEPENDENT COMPUTATION BEGINS HERE 
C 

20 IF(NN.EQ.l .AND. J.EQ.JAICALL OUTPUT 
K = K*1 
J = J*1 

IF(K.LT.IOO) GO TO 40 
0T=STAB*10. 

IE(K.LT.200) GO TO 40 
DT = 100. 

DO 30 N = 1,NA 

0T3 = (R(N*l)-R(Nn**2/ALPHA/4. 

IF (D1 .LT.DT3) GO TO 30 
DT = DT3 
30 CONTINUE 
40 TIME = TIME*OT 
UN(1) = 0. 

C 

C IN-'ERIOR POINT COMPUTATION BEGINS HERE 

C 

LOOP = 0 

50 DO 70 N = 2.NA 
NO = N»LOOP 
NM = NQ-^1 

UX = (U(NO)-U(NM) )*DDX 

UXX = (U(N*1)*U(N-1)-2.*U(N))*DDX2 

AE = AB(N)»STS*AC(N)/S2 

AG = AO(N)/S2 

UT = AE«UX+AG*UXX 

IF (LOOP.EO.l) GO TO 60 

UN(N) = U(N)»UT"DT 

GO TO 70 

60 UN(N) = (UO(N)+U(N)*UT«DT)*.5 
70 CONTINUE 
C 

C BOUNDARY POINT COMPUTATION BEGINS HERE 

C 

OTl = DT 
TIM = 0. 

TI = TIME-DT 
KIM=30 

OTl = DT/FL0AT(K1M) 

80 CONTINUE 

IF (LOOP.EQ.O) SN = S>ST*DT 
IF (LOOP.EO.l) SN = .5*(S*S0*ST*DT) 

TNAMN = UN(NA-1)/SN/Z(NA-1) 

IF (LOOP. EQ. 01 GO TO 90 
T(NC) = UO(NC)/SO 
U(NC) = UO(NC) 

S = SO 

90 TI=TIME-DT/2. 

P=P0*(1 .♦.859*DPC*AJ1*C0S(P0MEGA*TI-THETA) ♦.859*0PCl*AJl*C0S(2.*Pf 
1MEGA*TI) ) 

TF = TF0*(P/P0)**GA 
TAV = (TF*T(NC) )*.5 
DO = D(P»TAV) 
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• CO n ?.(j 
001123 
001125 
001132 
001140 

001161 

001200 

001204 

001206 

001224 

001230 

001237 

001244 

001251 

001257 

001261 

001262 

001266 

001271 

001274 

001277 

001300 

001301 

001305 

001307 

001310 

001325 

001332 

001336 

001341 


rOAV = P*EMB/ (RR*TAV) 

VOR = VIS/ROAV 

DDMl = .6* (VOR/DD)**. 333333 

AS = SORT (GAMma*32.2*P/ROAV) 

VR=CC1»AS* (COS (P0MEGA*TI«PID2-THETA) ♦0PC1/DPC*C0S(2.*P0MEGA*TI*PI1 
12 ) ) 

VH=CC2»AS* (COS (P0MEGA»TI-THETA) ♦0PC1/DPC*C0S(2.*P0MEGA*TI) ) 

VGA = VGAF* (1 ,-(S/SIT)o«3) 

OV = VGA-VD 

VT =: 5.65*(R0AV*ABS(0V) )*o.l6*(VIS/S)«*.84*DV/RH0L/S 
VO=VD*VT*DT/2. 

VEL--SORT(OV**2*VR**2*VH**2) 

DOMl = SORT (2.»VEL*S/VOR) 

H = AKR/2./S*(2.*DUM2*OOMl) 

AKG = (2.*0UM1*D0M1)*DD*0UM4/S/TAV 

AKGP=AKG*P 

S1=S 

95 PPV=PV (T (NO ) 

DUM=1 ,-PPV/P 

IE (0UM.6T.1 .E-20) GO TO 100 
WRITE (6»310) 

CALL OUTPUT 
CALL EXIT 

100 DUM=AKGP*ALOG(OUM) 

STN = OUM/RHOL 
DUM = -DUM 

TR = -( (T(NC)-Tn*H-(T(NC)-TFO)*DUM*CPV*OUM*SL(T(NC) ) >/AK 

TNAM = T(NA-1 j*(TNAMN-T(NA-in«TIM/nT 

TX=(T(NC)-TNAM)/TW00X 

SIM =Sl*ST«OTl 

RNA= S1*Z(NA) 


UU I 

001345 

001351 

001355 

001363 

001367 

001371 

001372 

001373 

001375 

001377 

001402 

001404 

001406 

001410 

001414 

001420 

001422 

001423 

001424 

001425 

001426 

001430 

001431 

001432 


no 


120 


C 

C 


IPNA = TX»YNA/S1 

TRR = (TR-TRNA)/(S1-RNA) 

UR = T(NC)+S1*TR 

UT = ALPHA* (S1*TRR*2,*TR) ♦UR*ST 

UN(NC) = U(NC)*UT*DT1 

T(NC) = UN(NC)/SIM 

ST = STN 

Sl= SIM 

TIM = TIM*DT1 

IF (TIM.LT.DT)GO to 95 

IF (KIM.GT.O) SN = SIM 

no 110 N = 1,NC 

UO(N) = U(N) 

R(N) = SN*Z(N) 

IF (N.NF.n T(N1 = UN(N)/R(N) 

U(N) = UN(N) 

SO = S 
S = SN 
S2 = S**2 
ST = STN 
STS = ST/S 

IF (LOOP.EQ.l) GO TO 120 

LOOP =1 ■ 

GO TO 50 
CONTINUE 

INTERIOR POINT CALCULATION AND BOUNDARY POINT CALCULATION END HERl 
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001432 

001435 

001436 

001440 

001441 

001442 

001444 

001445 

001445 

001456 

001460 

001463 

001465 

001467 

001477 

001500 

001505 

001510 

001512 

001515 


001515 

001523 

001527 

001531 

001532 

001533 

001536 

1/ O A 

O0lb3<* 

001560 

001565 

001567 

001574 

001577 

001577 

001601 

001602 

001607 

001610 

001616 

001621 

001622 

001622 

001624 

001624 

001624 

001630 

001632 

001633 

001642 

001643 

001646 

001653 

001654 


C 

wnoT = OUr^o.«S?*5T 
KKl = KKl+1 

IF (KK1.LT.3) GO TO 130 
KK2 = KK2+1 
KKl = 0 

WWW(KK2) = -WnOT 
TTT(KK2) = TIME 
130 CONTINUE 

T(1 )=(UN(2)eZ(3)**2/7(2)-UN(3)«Z (2) **2/Z ( 3) ) / (Z ( 3) **2-Z ( 2) **2) /SN 

IF (NN.EO.l) GO TO 140 

PATIO = (S/SIT)**3 

IF (RATIO .GE. 0.1 ) GO TO 20 

TR=TIME-TINT 

WRITE (6»350) NNfRATIO 

GO TO 150 

140 IF (S.GT.0.4641589»SIT) GO TO 20 
WRITE (6»360) 

TB = TIME 

IF (NN.EO.l) T90 = TB 
150 CONTINUE 
C 

C COMPUTATION OF SUMMATION HISTORIES BEGINS HERE 

C 

DUMAl sTB+OTBAR* (NN-1 ) 

0UMA2=TAU+0TBAR*(NN~1 ) 

IF(TB.GE.TAU)60 TO 155 
MARK=1 

DO 153 KK=1»NY 
IF(NN.E0.1)W(KK)=0. 

TK=TB* (KK-1 )*0TST 

IF(TK.GT.OUMA1.AND.TK,LT.OUMA2>TK=-100. 

IF(TK.GT.T8+DUHA2) TK=-100. 

IF(TK.GE.DUMA2) TK=TK-TAU 
INDEX=KK2 

W(KK)=FW00T (TK) *W(KK) 

153 CONTINUE 
GO TO 220 

155 DO 210 KK=1»NY 
MARK=1 

TK = TB*(KK-1)*DTST 
1 = 1 

160 T1 = TB+ (I-l ) «OTBAP 

IF (TK.GT.Tl)’ GO TO 170 
II = I 
GO TO 180 
170 I = !♦! 

GO TO 160 
180 CONTINUE 

IF (NN.EO.l) W(KK) = 0. 

INDEX = KK2 
DO 190 N = 1*2 

IF (NN.LT.II.AND.N.EQ.l) GO TQ 190 

JJ = N-1 

TI = TK-JJ»TAU 

W(KK) = W(KK)*FWDOT(TI) 

MARK=MARK»1 
190 CONTINUE 
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001656 

001661 

001663 

001671 

00167A 

001677 

001704 

001705 

001710 

001712 

001715 

001720 

001733 

001734 

001736 

001746 

001747 


001755 

001756 

001757 

001760 

001761 

001762 

001763 

001764 

001765 

001766 

001767 

001770 

001771 

001772 

001777 

002001 

002011 

002030 

002040 

002050 

002067 

002077 

002107 

002112 

002115 

002120 

002123 

002126 

002131 

002134 

002137 

002143 

002147 

002153 

002157 

002166 

002167 

002171 


NT = TK/TAU 

DO 200 NT = 2. NT 

T2 = NI*TAU* (NN-n *0TB4R 

ir (TK.LT.T2) GO TO 200 

TI = TK-NI«TAIJ 

W(KK) = W(KK)*FWOOT<TI) 

MARK=MARK*1 
200 CONTINUE 
210 CONTINUE 
220 CONTINUE 

WRITE (6.380) 

WRITE (6,370) (W(KK).KK = l.NY) 

WAV = 0. 

DO 230 KK = l.NX 

230 WAV = WAV+(W(KK)*W(KK*1))/2.*0TST 
WAV = WAV/TAU 

WRITE (6.390) WAV 
C 

C COMPUTATION OF RESPONSE FACTORS BEGINS HERE 

C 

PSQ=0. 

P1SQ=0. 

P2SQ=0. 

PISO=0. 

PI1SO=0. 

PI2SQ=0. 

RNL=0. 

R1=0. 

R2=0. . 

AINL=0. 

AH-0. 

AI2=0. 

00 231 KK=1.NY 
TK=TB* (KK-1 )*0TST 

1 ♦ 

IF(KK.EO.l.OR.KK.EO.NY) PM=.5 

APSQ=.859*AJ1*(0PC*C0S(P0MEGA«TK-THETA) ♦DPC1*C0S(2.*P0MEGA*TK) ) 

AP1S0=.859*AJ1*0PC*C0S(P0MEGA*TK-THETA) 

AP2S0=.859»AJ1*DPC1*C0S(2.*P0MEGA*TK) 

APIS0=,859»AJ1 *(0PC*SIN(P0MEGA*TK-THETA) ♦0PC1*SIN(2.*P0MEGA*TK) ) 

AP11SQ=.859*AJ1*0PC*SIN(P0ME6A«TK-THETA) 

AP12S0=,859*AJ1*0PC1*SIN(2.»P0MEGA«TK) 

WPRM=(W(KK)-WAV)/WAV 

RNL=RNL*PM*WPRM*APSQ 

R1=R1+PM*WPRM*AP1S0 

R2=R2*PM»WPRM«AP2SQ 

ainl=ainl+pm«wprm*apisq 

All=All*PMOWPRM»APIlSQ 

AI2=AI2*PMoWPRM*API2SQ 

PS0=PS0+PM*APS0**2 

P1SQ=P]SO*PM*AP1SO**2 

P2SQ=P2SQ+PM*AP2S0**2 

PISQ=PISQ+PM*APISQ**2 

PIlSQ=PIlS0+PM*APrlSQ»*2 

231 PI2S0=PI2S0*PM*API2S0**2 
RNL=RNL/PSQ 
R1=R1/P1S0 
R2=R2/P2SQ 
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00217? 

00217A 

002175 

002177 

002225 

002243 


002270 

002273 

002275 

002277 

002300 

002301 

002302 

002335 

002336 

002337 

002340 

002341 

002344 

002352 

002354 

002356 

002377 

002413 

002420 

w w C T»_ * 

002425 

002426 

002426 

002426 


002426 


002426 


002426 

002426 

002426 

002426 


002426 


ainl=ainl/piso 

AI1=AI1/PI1SQ 

AI2=AI2/PI2SQ 

C APF=0MEGA*(STT/. 000 166) **1.5* ( 800. /VGAF)««.3333*(.1/0PC) “*.3333* 

1 (43200. /PO)**. 3333 

CAPG=CAPF» ( VGAF/800. ) ** (7./3. ) * (DPC/. 1 ) **.633 
WRITE < 6,400 )RNL» R1 ,R2,AINL,AI1 , AI2,CAPF,CAPG,T90 
C 

C CALCULATION OF N AND TAU BEGINS HERE 

C 

AP=.859*AJ1*DPC 

8P=.859*AJ1*0PC1 

WT=PI2/2000. 

SAVEX=0. 

SAVJ=0. 

DO 240 1=1,2000 

EX=RNL/AINL* (AP**2*SIN(WT) ♦BP**2*SIN(2.*WT) ) *BP**2* ( 1 .-COS ( 2. *WT ) ) 
1»AP**2*(1 .-COS(WT) ) 
lE(FX) 232,238,233 

232 J=-l 

GO TO 234 

233 J=1 

234 IF( J*SAVJ)239,238,239 

238 WTS=WT+P 12/2000.* (EX/ (SAVEX-EX) ) 

PTAU=WTS/POMEGA 

ANGLE=WTS*360./PI2 

PEN=RNL* ( AP**2*8P**2) / ( AP**2* (1 .-COS (WT) > ♦BP**2* ( 1 .-COS (2.*WT) > ) 
WRITE(6,500) WTS, ANGLE, PTAU, PEN 

239 WT=(I*1 )/2000.*PI2 
SAVFX=FX 

L.'TKI JmVJ = J 

CALL EXIT 
C 

250 EORMAT (1415) 

260 FORMAT (8E10.4) 

270 FORMAT (IHO, lOHRUN NUMBERI5,5X,2H0NI3, 1H/I2 , 1 H/I2// 

1 13H OUTPUT EVERYI6,6H STEPS//14H GA£ 

2 PRESSURE=E13.5,7H LB/FT2/17H GAS TEMPERATURE=E13.5,8H PANKINE/24F 

3 DROPLET INITIAL RADIUS=E13.5,3H FT/21H DROPLET TEMPERATURE=E13.5, 
48H RANKINE/llH FINAL VGA=E13.5,7H FT/SEC/18H DROPLET VEL0CITY=E13. 
55, 7H FT/SEC/51H AUTOMATIC STOP WHEN 90 PERCENT MASS HAS EVAPORATED 
6) 

280 FORMAT (9H0 ALPHA, 7X ,4HRH0L,6X, 12HGAS CONST ANT , 3X , 3HPC A , 9X , 3HPCF 
1 ,9X,3HTCA,9X,3HTCB/7E12,4//4X,2HMA,16x,2HMB,10X,2HMV,11X,1HA,11X,K 
2R,9X,3HSCH,9X,2HPP/7E12.4//4X,2HKB,8X.5H(CP) V,7X,5HKAPPA,7X,5HGAMF 
3A,7X,3HVIS/5E12.4/) 

290 FORMAT (35H PV (T)=EXP ( 16.928-1476.5/ (T-3. 568) ) /36H LAMBDA (T ) =61 . 3:* 
1+.5916«T-.00248*T**2/83H D (P,T) = (PCA*PCB) **1 /3/SQRT (MA*MB/ (M A *MB ) ) 
2*(TCA*TCB) **5/12/P*A* (T/SORT (TCA*TCB) )**B) 

300 FORMAT (36H0NUMBER OF INTERVALS INSIDE DR0PLET=I3,23H, STRETCHING 
1PARAMETER=E12.4/21H STABILITY PARAMETER=E12.A) 

310 EORMAT (6H0EAIL1) 

320 EORMAT (IHl, 36HPR0GRAM DROPS AUGUST 1972 VERSION ) 

330 FORMAT (46H THE NUMBER OF CALCULATED DROPS PER PERIOD IS I5/58H THE 
1 NUMBER OF ARTIFICIAL DROPS BETWEEN EACH REAL DROP IS I5/51H THE E 
2VAP0RATI0N FROM THE ENTIRE ARRAY IS SUMMED AT,I4,18H POINTS PER PE 
3RI0D) 

340 FORMAT(9HO OMEGA, 7X,5HTHETA,5X,10HDPC PK-PK,6X,3HAJ0,9X.3HAJ1 , 
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oo?':?.e, 

002426 

002426 

002426 

002426 

002426 


002426 

002426 

002426 


19X,3HAJ2»9X»3HRWR/7E12.4) 

350 FORMAT(1HO*5HOROP I5«36H CALCULATION ENDS WITH A MASS RATIO=E10.3>' 
360 format (1H0.5X02H90 PERCENT EVAPORATION COMPLETED) 

370 FORMAT (1X,10E10.3) 

380 FORMAT (1H1»6H W IS> 

390 F0RMAT(1X«5H WAV= E12.5) 

400 F0RMAT(1X»5H RNL= E12.5/1X»5H Rl= E12.5/1X»5H R2= E12.5/1X»5H INL= 
1 E12.5/1X,5H 11= E12.5/1X»5H 12= E12.5/1X»5H F = E12.5/1X.5H G = 
2E12.5/1X»5H T90=E12.5» 8H SECONDS) 

410 F0RMAT(1X»5H0PC1=E12.5»5X,9H(DP1/DPC) ) 

500 FORMAT (/// 9H WTS = E12.5*8H RADIANS/9H ANGLE = E12.5»8H DEGREES 

1/9H TAU = E12.5»8H SEC0NDS/9H N = E12.5) 

END 


1 



000003 


FUNCTION FWnOT (T) 

COMMON WWWOOOO) »TTT (3000) .W(400> * INDEX ♦ KK2 * VEL ♦ VO» AKG. DO » PPV t Hf 00 
IMl ,STT*S»TIME, J«K« JA,NA«NC»DX,OT*ST*SN»STN»WDOT.P,TF 
000003 COMMON SAVW (30*400) »MARK,KK»NAPT»NN.NP 

C 

000003 IF(T.LT.TTT(1) )FWDOT=0. 

000006 IF(T.LT.TTTd)) GO TO 30 

000010 IF (T .GT.TTT (KK2) ) FWDOT = WWW(KK2) 

000015 IF(T,GT.TTT(KK2) ) GO TO 30 

000021 DO 10 I = 1*INDEX 

000022 L = TNDEX+l-I 

000024 IF (TTT(L).GE.T ) GO TO 10 

000027 LM = L 

000030 GO TO 20 

000031 10 CONTINUE 

000034 20 LP = LM+1 

000036 INDEX = LP 

000037 EPS = (TTT(LM)-T ) / (TTT (LM) -TTT (LP) ) 

000044 FWDOT = WWW (LM ) ♦EPS* ( WWW (LP) -WWW (LM) ) 

C 

C BEGINNING OF ARTIFICAL DROP CALCULATION 

C 

000051 30 IF(NART.EO.O) RETURN 

000054 WW=FWOOT 

000056 OLAST=0. 

000057 IF(NN.NE.l) GO TO 40 

000061 SAVW(MARK*KK)=WW 

000065 ART=0. 

C 

V. «uL)S IN CONTRIBUTION FROM DROP 1 FOR INTERPCLATIGI, mIT,; LAST LAS." 

C 

000065 DO 35 I=1*NART 

000067 35 ART=ART+FLOAT(I)/FLOAT(NART^l)*WW 

000100 FWDOT=FWDOT+ART 

000101 RETURN 

000102 40 IF(NN,EQ,NP) DLAST=1, 

000106 ART=0. 

000107 DO 50 I=1«NART 

000111 50 ART=ART+SAVW(MARK*KK)^FL0AT(I)/FL0AT (NART+1)*(WW-SAVW(MARK»KK) ) 

1 ♦OLAST* ( WW-FLOAT ( I ) /FLOAT (NART+ 1 ) *WW) 

000136 SAVW(MARK*KK)=FWDOT 

000142 FWDOT=FWOOT*ArT 

000143 RETURN 

000144 END 
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00000 ? 

000002 

000002 


000002 

000004 

000021 

000024 

000025 

000027 

000031 

000053 

000054 

000055 

000060 

000073 

OOOlll 

000121 

000121 

000122 

000123 

000123 

000123 

000123 

O00123 

000123 


SUBROUTINE OUTPUT 

COMMON WWW (30001 »TTT (30001 »VU4001 ♦!N0EX»KK2»VEL»V0»AKG»00»PPV»H.n0 
lMl»SIT,S»TIME.J»K»JA»NA«NC.DX,DT»ST.SN,STN»WDOT*P»Tr 
COMMON SAVW(30«400) *MARK»KK*NART*NN*NP 

COMMON X(64) »Z(64) »AB(641 ♦AC(64) ♦AD(64) ,U0(64) .U(64) ,UN(64) ,T(64) » 
1R(64> 

OUM = TIME 

WRITE (6v40l K«0UM»SN»ST»W00T 
NH = <NC*ll/2 
NM = NC/2 
DO 10 N = 1«NH 
L = NM*N 

10 WRITE (6«30) N«R<N)»T(N)»L«R(L)*T(L) 

IE (K.EO.O) GO TO 20 
REY = 00M1**2 
DUM = (S/SIT)**3 
WRITE (6*501 H«REY*0UM»DT 
WRITE(6»70) VEL*VD*AKG«DO*PPV 
WRITE (6*60) P*TF 
20 CONTINUE 
J = 0 
RETURN 

30 FORMAT ( 19 *2E1 5.5 (I9*2E15.5) 

40 FORMAT (1HO//5X*4HSTEPI6*10X*5HT1ME=E12.4*5X»3HRS=E12.4/24X,6HI?SOO 
1T=E12.4*5X*5HW00T=E12.4//18X*1HR.13X.1HT*24X*1HR*13X,1HT) 

50 format (2X»2HH=E10.3*1X*4HREY=E10.3*1X*11HMASS RATIO=E10.3*1X«3HOT 
1=E10.31 

60 format (2X*2Hp=E10.3*1X*3HTF=E10.3) o 

70 FORMAT (2X*4Hy£L=E 10. 3*3HVO=E10,3*4HAKG=E1 0,3* 3HUl)=tl0.3*.‘*MHKV=ti0. 
13) 

END 
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ro = 50 MICRONS 



fj. = 300psi 
AP, =0.2 Ft 
AI|=0.8AF| 

f =3000 CPS 
V^j=IOOFT./SEC. 

Ugf = l200FT./SEC. 


TIME (SEC.) 



p 




FIG. 2 ARRAY EVAPORATION RATE VS. TIME 
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10 ^ 2 4 6 8 10 * 2 

f- FRE( 

FIG. 3 RESPONSE FACTOR VS. FREC 
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Pc = 300 psia 
Vh = I00FT7SEC 
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f - FREQUENCY 

FIG. 5 RESPONSE FACTOR VS. FREQUENCY: EFFECT OF FINAL GAS 
VELOCITY 







AP, =0.2 Pc 
Ap2=0.8AP, 
Vn=l00 FT SEC. 
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FINAL GAS VELOCITY 
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'*wf(Tr(-5^r 

FI6.7 IN PHASE NON-LINEAR RESPONSE FACTOR VS. FREQUENCY FACTOR 
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FIG. 8 IN PHASE FUNDAMENTAL RESPONSE FACTOR VS. FREQUENCY 
FACTOR 
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FI6.9 FIRST HARMONIC RESPONSE FACTOR VS. FREQUENCY FACTOR 
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FREQUENCY FACTOR 
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PHASE RESPONSE FACTOR VS. 
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FIG. 12 RRST HARMONIC OUT OF PHASE RESPONSE FACTOR VS. 
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